Exact Criterion for Determining Clustering vs. Reentrant Melting Behavior for 

Bounded Interaction Potentials 
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We examine in full generality the phase behavior of systems whose constituent particles interact 
by means of potentials which do not diverge at the origin, are free of attractive parts and decay 
fast enough to zero as the interparticle separation r goes to infinity. By employing a mean field- 
density functional theory which is shown to become exact at high temperatures and/or densities, 
we establish a criterion which determines whether a given system will freeze at all temperatures or 
it will display reentrant melting and an upper freezing temperature. 
PACS: 61.20.-p, 61.20.Gy, 64.70.-p 
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I. INTRODUCTION 



The phase behavior of systems whose constituent par- 
ticles interact by means of pair potentials diverging at 
the origin is a problem that has been extensively stud- 
ied in the last few decades. The whole range of inverse 
power-law pair potentials have been examined, ranging 
from hard spheres (HS) to the one-component plasma 
(OCP) and it has been established that excluded vol- 
ume effects are mainly responsible for bringing about the 
freezing transition. The crystal structure in which a liq- 
uid freezes is subsequently determined by the steepness of 
the repulsion, with hard repulsions favoring a face cen- 
tered cubic (fee) lattice and soft ones a body centered 
cubic (bec) lattice |jj . Power-law diverging potentials re- 
sult into freezing at arbitrarily high temperatures. How- 
ever, the divergence of the potential alone is not enough 
to cause such a phenomenon, as demonstrated recently 
by Watzlawek et al. who employed a logarithmically 
divergent pair potential, suitable to describe effective in- 
teractions between star polymers in good solvents. It 
was shown that the strength (prefactor) of the logarith- 
mic potential, determined by the number of arms / of the 
stars, is crucial in determining whether the system freezes 
or remains fluid at all densities. As a consequence, the 
phase diagram of star polymers was predicted to display 
reentrant melting and a critical freezing value f c = 34 of 
arms, such that for f < f c the system remains always 
fluid @. 

Another interesting class of interactions are those 
which do not diverge at the origin, i.e., they are bounded. 
Such potentials arise naturally as effective interactions 
between the centers of mass of soft, flexible macro- 
molecules such as polymer chains ||, dendrimers j|, 
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polyelectrolytes etc. Indeed, the centers of mass of 
two macromolecules can coincide without violation of 
the excluded volume conditions, hence bringing about 
a bounded interaction. Moreover, the same mechanisms 
that exist for tuning the usual, diverging interactions be- 
tween colloidal particles can be applied in order to tune 
the bounded interactions: the solvent quality, tempera- 
ture, chain length, salt concentration etc. will all affect 
the effective potential. Thus, it appears to be useful to 
consider such potentials in some generality in order to 
be able to draw conclusions about the expected phase 
behavior of systems interacting by means of these. 

Two model systems in this category have already been 
studied in some detail. One is the penetrable spheres 
model (PSM) || r|, in which the pair potential is a pos- 
itive constant e for distances r < a and vanishes oth- 
erwise. The other is the Gaussian core model (GCM), 
introduced in the mid-seventies by Stillinger J/J. In 
the GCM, the pair potential v(r) has the form v(r) = 
e exp[— (r/<r) 2 ], with e being an energy and a a length 
scale. It has been shown that the GCM models very ac- 
curately the effective interactions between the centers of 
mass of linear polymer chains §, § |, [k], [ll], |§, p[ @. 

The PSM was studied by means of cell-model calcula- 
tions and computer simulations liquid-state integral 
equation theories |l5| and density- functional theory [jl6| ; 
the fluid structure of the PSM has been further stud- 
ied recently by Rosenfeld et al. by using ideas based 
on the universality of the bridge functional [[l7| . It was 
found that no reentrant melting takes place because the 
solid always lowers its free energy by allowing for mul- 
tiply occupied crystal sites, a mechanism that is called 
clustering Q . The clustering mechanism stabilizes there- 
fore the solid at all temperatures. Hence, the topology 
of the phase diagram of the PSM is similar to that of 
power-law diverging potentials, when details about the 
clustering structure of the solids are disregarded. On 
the other hand, the GCM has been studied even more 
extensively by means of molecular dynamics simulations 
fll8|, |l9[ , high-temperature expansions pfj and the discov- 
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ery of exact duality relations in the crystalline state . 
Recently, a full statistical-mechanical study of the GCM 
was performed and it was established that the topology 
of the phase diagram of the GCM resembles that of star 
polymers. Freezing and reentrant melting accompanied 
by an upper freezing temperature were quantitatively cal- 
The question that arises, therefore, is the 



culated 

following. Given a nonattractive and bounded pair po- 
tential which satisfies the following requirements guaran- 
teeing stability and the existence of the thermodynamic 
limit @: 

(i) it is bounded; 

(ii) it is positive definite; 

(iii) it decays fast enough to zero at large separations, 
so that it is integrable and its Fourier transform 
exists; 

to which topology belongs the phase diagram of the sys- 
tem? In this paper, we present an exact criterion which 
gives an answer to this question and show representative 
results for model systems which confirm its validity. The 
rest of the paper is organized as follows: in section II 
we present the physical arguments supporting the mean- 
field theory of the models and in section III we discuss 
the existence of a spinodal instability in this theory and 
its implications on the phase bahavior. We present a 
systematic comparison between theory and simulation in 
section IV and we draw the generic phase diagrams of 
such systems in section V. Finally, in section VI we sum- 
marize and conclude. 



II. THE MODEL AND THE MEAN-FIELD 
LIMIT 

We will work with a general interaction v(r) — s(f>(r/a) 
satisfying the requirements put forward above. Here, e 
and a are an energy and a length, respectively, and 4>{x) 
is some dimensionless function. The latter does not have 
to be analytic, i.e., discontinuities in the potential or its 
derivatives are allowed. Without loss of generality, we 
assume 0(0) = 1. Let us call (j)(Q) = a~ 3 (j){Q) the di- 
mensionless Fourier transform (FT) of the interaction. 
For more concreteness (and for the purposes of demon- 
stration) we introduce in addition the family of bounded 
potentials (r) depending on a tunable parameter £, 



1 + e~ a ^ 



Vf(r) = e -, 



(1) 



where £ is a 'smoothing parameter' having dimensions 
of length. The case £ = recovers the PSM whereas 
as £ grows the interaction becomes smoother. Due to 
its resemblance to the Fermi-Dirac distribution, we call 
this family the Fermi distribution model (FDM). The 
additional factor 1 + e~ a '^ in the numerator of the rhs of 



cq. (pj) ensures that the potential varies from e at r = 
to zero at r — > oo, for all £. 

We introduce dimensionless measures of temperature 
and density as 



t 
V 



6 P<T = ^ 



(2) 
(3) 



where fcs is Boltzmann's constant and p — N/V is the 
density of a system of N particles in the volume V. We 
will refer to 77 as the 'packing fraction' of the system. 

The key idea for examining the high temperature 
and/or high density limit of such model systems is the 
following. We consider in general a spatially modulated 
density profile p(r) which does not vary too rapidly on 
the scale a set by the interaction. At high densities, 
pa 3 3> 1, the average interparticle distance a = p^ 1 / 3 
becomes vanishingly small, and it holds o < <r, i.e., the 
potential is extremely long-range. Every particle is si- 
multaneously interacting with an enormous number of 
neighboring molecules and in the absence of short-range 
excluded volume interactions the excess free energy of the 
system can be approximated by a simple mean-field 
term, equal to the internal energy of the system: 



drdr'«(|r-r'|)/)(r)p(r / ), (4) 



with the approximation becoming more accurate with in- 
creasing density. Then, eq. (^) immediately implies that 
in this limit the direct correlation function c(r) of the 
system, defined as p3] 



c(|r-r'|;p) 



5 2 pF cx [p(r)} 
p(t)^p 8p(r)5p(r') ' 



lim 



(5) 



becomes independent of the density and is simply pro- 
portional to the interaction, namely 



c(r) = —(3v(r). 



(6) 



Using the last equation, together with the Ornstein- 
Zernike relation p5j , we readily obtain an analytic ex- 
pression for the structure factor S(Q) of the system as 



S(Q) 



1 



l + pt-^iQ) 



(7) 



This mean-field approximation (MFA) has been put for- 
ward and examined in detail in the context of the Gaus- 
sian core model independently by Lang et al. p2| and by 
Louis et al. ]2l| . The model is particularly relevant from 
the physics point of view, due to its connection to the 
theory of effective interactions between polymer chains 
§, |j. Here, we establish the validity of the MFA at 
high densities for bounded, positive-definite interactions 
in general and we examine its implications for the global 
phase behavior of such systems. 
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Bounded and positive-definite interactions have been 
studied in the late seventies by Grewe and Klein [^7], p8| . 
The authors considered a slightly different model than 
the one considered here, namely a Kac potential of the 
form: 



v(r) = j tpi'yr), 



(8) 



where d is the dimension of the space and 7 > is a 
parameter controlling the range and strength of the po- 
tential. Moreover, ip( x ) is a nonnegative, bounded and 
integrable function: 

< rp(x) < A < 00, C = / d d xij;(x) < 00. (9) 



Grewe and Klein were able to show rigorously that at the 
limit 7 — > 0, the direct correlation function of a system 
interacting by means of the potential (Q) is given by eq. 
(|J) above. The connection with the case we are discussing 
here is straightforward: as there are no hard cores in the 
system or a lattice constant to impose a length scale, the 
only relevant length is set by the density and is equal to 
p-!/3 m our mo del and by the parameter r y^ 1 in model 
(H). In this respect, the limit 7 — ► in the Kac model 
of Grewe and Klein is equivalent to the limit p — » 00 
considered here. However, in the Kac model the strength 
of the interaction goes to zero simultaneously with the 
increase in its range. Moreover, the validity of the mean- 
field expression at large but finite densities and at low 
temperatures has not been tested in detail. 



III. SPINODAL INSTABILITY AND FREEZING 

We employ the MFA as a physically motivated work- 
ing hypothesis for now and, by direct comparison with 
simulation results, we will show later that it is indeed 
valid. Within the framework of this theory, an exact cri- 
terion can be formulated, concerning the stability of the 
liquid phase at high temperatures and densities. The 
function <f>(x) was assumed to be decaying monotonically 
from unity at x = to zero at x — > 00. For the function 
4>(Q), there are two possibilities: (i) It has a monotonic 
decay from the value <j>(Q = 0) = cr -3 J dx(f>(x) > to 
the value 4>{Q) — at Q — > 00. We call such potentials 
Q + -potentials. Obviously, the Gaussian interaction be- 
longs to this class, (ii) It has oscillatory behavior at large 
Q, with the implication that it is a nonmonotonic func- 
tion of Q, attaining necessarily negative values for certain 
ranges of the wavenumber. We call such potentials Q - 
potentials. Long-range oscillations in Q-space imply that 
<fi(x) changes more rapidly from unity at r — to zero at 
r — > 00 in the Q ± -class than in the Q + one. Moreover, 
let us call Q* the value of Q at which <j>{Q) attains its 
minimum, negative value. 

If we are dealing with a (^-potential eq. (0) implies 
that S(Q) has a maximum at precisely the wavevector 
Q» where 4>{Q) attains its negative minimum, — |<^(Q*)| 



and this maximum becomes a singularity at the 'spinodal 
line' pt~ 1 \4>(Q tf )\ — 1, signaling the so-called Kirkwood 
instability of the system |2^, |2^, EnJ. The theory 
has a divergence, implying that the underlying assump- 
tion of a uniform liquid is not valid and the system must 
reach a crystalline state. Indeed, on the basis of the 
fluctuation-dissipation theorem, S{Q) can be interpreted 
as a response function of the density to an infinitesimal 
external modulating field at wavenumber Q [£5| and a di- 
verging value of this response function clearly signals an 
instability. If the Fourier transform of <j)(x) has negative 
Fourier components, then an increase in temperature can 
be compensated by an increase in density in the denom- 
inator of eq. (0), so that S(Q*) will have a divergence 
at all t. We thus conclude that Q^-systems freeze at all 
temperatures. 

If we are dealing with a Q + -potential (4>(Q) mono- 
tonic), then eq. (Q) implies that S(Q) is also a monotonic 
function of Q at high densities |2^| . For such potentials, 
one can always find a temperature high enough, so that 
the assumptions of eq. (||) hold and then eq. @) forces 
the conclusion that freezing of the system is impossible 
at such temperatures. This does not imply, of course, 
that such systems do not freeze at all; one simply has 
to go to a low enough temperature and density, so that 
the mean-field assumption does not hold and the inter- 
action is much larger that the thermal energy. Then, the 
system will display a hard-sphere type of freezing, to be 
discussed more explicitly below. An upper freezing tem- 
perature t n must exist for <5 + -potentials, implying that 
such systems must remelt at t < t u upon increase of the 
density. Hence, we reach the conclusion that Q + -systems 
display an upper freezing temperature and reentrant melt- 
ing. The criterion says nothing about the crystal struc- 
ture of the solid, however, which always depends on the 
details of the interaction as well as the density |2^, |3^] . 

For potentials is in the <5 + -class, the mean-field ar- 
guments presented above hold not only at high temper- 
atures but also at low ones, provided that the require- 
ment pa 3 ^ 1 is satisfied, because these are molten at 
high densities for all nonzero temperatures. The valid- 
ity of the mean-field theory for Q + -type systems, even 
at very low temperatures, was confirmed recently by di- 
rect comparison with simulation results for the particular 
case of the Gaussian potential (22) . If the potential is in 
the Q ± -class, the mean field approximation holds pro- 
vided that the system is not already frozen, as we will 
confirm shortly. Moreover, both kinds of systems display 
an unusual kind of 'high density ideal gas' limit. Indeed, 
taking the expression (0) for S(Q) and using the relation 
S(Q) = 1 + ph(Q) p5| , where h(Q) is the dimensionless 
Fourier transform of the pair correlation function h(r) of 
the uniform fluid, we obtain: 



HQ) 



t-^(Q) 



(10) 



l+pt-^Q) 

At low Q's, where <f>(Q) is of order unity, the term 
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proportional to the density in the denominator domi- 
nates in the limit of high densities and h(Q) scales as 
— 1/p — ► 0. At high Q's, the Fourier transform <f>{Q) in 
the numerator is itself small, with the result that h(Q), 
and hence also the correlation function h(r), is approach- 
ing zero. This, in turn, means that the radial distri- 
bution function g(r) — h(r) + 1 is very close to unity 
in this limit and deprived of any significant structure 
for all values of r and it only has some small structure 
at small r, which is in fact more pronounced for Q^- 
potentials than for Q + ones. In this limit, the hyper- 
netted chain (HNC) closure becomes exact, as the exact 
relation g(f) = exp[/?u(r) + h(r) — c(r) — B(r)], combined 
with the limits g(r) —> 1, h(r) — > and c(r) — > —f3v(r) 
forces the bridge function B(r) to vanish. Moreover, eqs. 
(0) and ([l0]) reveal that the systems obey a scaling law, 
namely that the functions S(Q) and th(r) do not depend 
on p and t separately but only on the ratio p/t. 

Systems in the Q ± -class freeze before the spinodal 
is reached. In order to make quantitative predictions, 
we invoke the empirical Hansen- Verlet freezing criterion 
p3| , p4| , which states that a system crystallizes when 
S(Q) at its main peak attains, approximately, the value 
S{Q*) = S m = 3. Although this criterion was originally 
put forward for hard, atomic interactions (HS, Lennard- 
Jones etc.), recent detailed analyses have demonstrated 
that it holds for the freezing and the remelting transi- 
tions of ultrasoft particles such as star polymers [p[ |35j| 
and even for the nondiverging Gaussian interaction p2| . 
Hence, we assume that it is valid for the general class of 
systems we consider here and combining it with eq. ([?]), 
we obtain the equation of the freezing line U (77) as 

*ffa) = y (Q g -L *? = 2 - 864 \kQ.)w (11) 

7T(1 - b m ) 

The value |0(Q*)| determines the slope of the freezing 
line at the high (t, 77) part of the phase diagram. 



IV. COMPARISON WITH SIMULATIONS 

We now wish to put these arguments in a strong test, 
using the concrete family of the FDM, given by eq. (Q). 
First of all, we have calculated the Fourier transform of 
the potential v^ir) of eq. ([!]) numerically, establishing 
that members of the FDM with £ < £ c belong to the 
class and members with £ > £ c = 0.49697 to the Q + one. 
The GCM is also a member of the latter class that we 
will discuss in what follows. 



A. Systems displaying clustering 

As examples of systems displaying clustering transi- 
tions, we have taken the extreme (and by now well- 
studied case) case £ = (the PSM) as well as the case 
£ = 0.1 of the Fermi distribution model of eq. ([!]). We 



have performed standard Monte Carlo (MC) NVT sim- 
ulations for a large number of values for the temperature 
and density. We begin with the PSM for which the ana- 
lytical expression (0) takes the form 



S(Q) 



1 + 24r)t~ 



sin(<3er) — (Qa) cos(Q<r) 



(12) 

The high temperature-high density freezing line of eq. 
( |TT| ) takes for this choice of £ the form tjjrj) — 1.033?7. 
To test the analytical expression of eq. (JI4), we move 
along the 'diagonal' t = 77, a combination that lies al- 
most on the Hansen- Verlet estimate for the location of 
the freezing line. In Fig. [l] we show the comparison of 
the analytical results with those obtained from the MC 
simulations for S(Q) and we also demonstrate that the 
MC curves for the quantity th(r) all collapse onto a single 
line, amply demonstrating the validity of the mean-field 
approximation for the PSM. In order to further inves- 
tigate the validity of the MFA, we have performed MC 
simulations in a variety of thermodynamic points and we 
present a selection of the obtained results. We present 
a selection of these in Figs. || and || and discuss them 
below. 

In Fig. ||(a) we show a comparison between the MC 
and MFA results for the radial distribution function 
g(r) along the 'diagonal' t = r). It can be seen that 
the agreement between the two is already very good at 
t = i] = 4.0 and thereafter it improves markedly with in- 
creasing temperature and density. The results obtained 
from the present theory are of the same quality as those 
obtained by Fernaud et al. [ fL5| , who used the sophisti- 
cated zero-separation (ZSEP) closure to investigate the 
liquid structure of the system. This closure involves three 
self-consistency parameters, determined in such a way 
that the virial-compressiblity, Gibbs-Duhem and zero- 
separation consistency conditions are fulfilled. At the 
same time, the present results are of the same quality 
as the recently obtained results of Rosenfeld et al. |[L7) , 
based on ideas of the universality of the bridge functional. 

In Fig. §(b) we perform the same comparison but now 
at fixed temperature t — 5.0 and increasing packing frac- 
tion r\. As can be seen, at this temperature, the MFA, 
which was originally formulated as a high-density approx- 
imation, proves to perform extremely well even at inter- 
mediate packings, 7/ = 0.5 for instance. This is a direct 
consequence of the boundedness of the interaction com- 
bined with a temperature t 3> e. Indeed, for small den- 
sities, the direct correlation function tends to the Mayer 
function, c(r) = exp [— f3v(r)] — 1 [^5|. If we are deal- 
ing with a bounded interaction at high temperature, we 
can linearize the exponential, obtaining c(r) = —(3v(r) at 
low densities, which matches with the MFA expression, 
eq. @, at high densities, thus leading to the conclusion 
that the MFA is an excellent approximation at all densi- 
ties. For unbounded interactions the linearization of the 
exponential is evidently impossible. 

In Fig. |3| we present a comparison between MC and 
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FIG. 1: (a) The product th{r) for a FDM with 
£ = (PSM), along the diagonal t = r\ at high 
packing fractions, as obtained from MC simulations. 
The results close to r — are noisy due to poor 
statistics there. All results collapse onto a single 
curve, (b) The corresponding structure factors S(Q), 
shown together with the analytical result of eq. (12). 



MFA at fixed packing fraction r\ = 3.0 and increasing 
temperature. As can be clearly seen, the validity of 
the MFA improves with increasing temperature. For 
bounded interactions, an increasing temperature implies 
a 'washing-out' of the correlation effects caused by the 
(increasingly weak) interaction effects and a tendency of 
the system towards the particular 'high-density ideal gas' 
limit characterized by the tendency of the function g(r) 
towards unity. However, it is an interesting peculiarity 
of these systems that unlike the usual ideal gas, the limit 
g(r) — > 1 (or, equivalently, h(r) — > 0) does not imply 
a corresponding limit S(Q) — > 1. Though the Fourier 
transform of h(r), h(Q), tends to zero as p , this is 
compensated by the large density p, so that the struc- 
ture factor S(Q) = 1 + ph(Q) displays the signature of 
strong ordering through the pronounced peaks seen in 
Figs. |l|(b) and|(b). 



FIG. 2: (a) The function g(r) of the PSM for 
selected points along the 'diagonal' t — rj as 
obtained from theory (thick lines) and simulation 
(thin lines). (b) Same but now for fixed tem- 
perature t — 5 and increasing packing fraction r/. 



Further, we performed MC simulations at selected 
points deeply inside the region t < tf (77) , finding that 
the obtained structure factors displayed Bragg peaks and 
hence confirming the prediction that the system is frozen 
there. Putting all our results together, we draw in Fig. 
^ a semi-quantitative phase diagram of the PSM, accom- 
panied by an assessment of the validity of the MFA at 
selected thermodynamic points. The MFA appears to 
be an excellent approximation at all densities above the 
temperature t = 3.0. Hence, we take as an estimate for 
the freezing line above t — 3.0 the MFA-Hansen-Verlet 
line tf = 1.03377 = 77; for lower temperatures, we sim- 
ply connect the point (77, t) = (3.0, 3.0) with the point 
(77, t) — (0.5,0), which obtains from the consideration 
that at t = the PSM reduces to the hard sphere system 
which is known to freeze at a fluid density 77ns — 0.5. The 
monotonic shape of the freezing curve for low tempera- 
tures arises from detailed considerations there, which can 
be found in Ref. 10. 
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FIG. 4: The phase diagram of the PSM, along with 
the points where the mean field theory brings excel- 
lent agreement with simulation (filled circles), fairly 
good agreement (empty squares) and no good agree- 
ment (x -symbols). These symbols should help delin- 
eate the domain of validity of the mean-field theory. 



First, let us consider the number of particles iV c in 
the fluid phase whose centers are, on average, within a 
distance a from a given particle. The number N c is given 
by the formula: 



FIG. 3: (a) Same as Fig. ||(b) but now for fixed 
rj = 3.0 and increasing temperature, (b) The struc- 
ture factors at the thermodynamic points of (a), com- 
parison between theory (lines) and simulation (points). 



Next we present in Fig. || a comparison for the FDM 
with £ = 0.1. For this choice of £, the Hansen- Verlet- 
based freezing line takes the form tf = 0.712?7. The se- 
lected points lie in the fluid region and the comparison 
indicates once more the excellent accuracy of the MFA 
both for g(r) and for S(Q). The radial distribution func- 
tion g(r) of this model is deprived of the jump at r = a 
seen in the PSM; the latter is caused by the disconti- 
nuity of the PSM potential there. However, a similarity 
between the g(r)'s of the £ = and £ = 0.1 models is that 
they both attain their maximum values at full overlaps 
between the particles r = and thereafter they decay 
rapidly, featuring a depletion region around r « a. This 
is a characteristic pointing to a strong clustering property 
in the fluid phase, a property thereafter inherited by the 
incipient thermodynamically stable crystal; the number 
of particles 'sitting on top of each other' and thereby oc- 
cupying the same crystal site scales linearly with density. 
In order to corroborate this claim, we can argue in two 
different ways, using the liquid as a reference point. 



N c = 1 + 4tt/9 / r 2 g(r)dr. (13) 
Jo 

In Fig. U we show the function ATrr 2 g(r) within a particle 
diameter a for a sequence of points along the freezing 
line of the PSM. As all these curves tend to a common 
limit with increasing density, the integral Air J Q r 2 g(r)dr 
tends to a constant and hence N c oc p at high densities, 
where the second term on the rhs of eq. ( |l3"| ) dominates. 

Second, we can use the wavevector at which the 
fluid structure factor has a maximum in order to estimate 
for the lattice constant a of the incipient crystal through 
the relation a oc Q^ 1 . For the models at hand, this 
maximum is entirely determined by the pair potential; 
unlike in usual fluids featuring diverging interactions, for 
which scales as p 1 / 3 , in our case knows nothing 
about the density. Thus, all post-freezing crystals have 
the same lattice constant, although their average density 
is a linear function of the temperature. This clearly shows 
that clustering must take place in the crystal: by allow- 
ing more and more particles to occupy the same lattice 
sites, a practically constant effective density of clusters 
is maintained in the crystal, thus leading to a density- 
independent lattice constant. 
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FIG. 5: (a) Comparison between theory (thick lines) 
and simulation (thin lines) results for g{r) of a FDM 
system with £ = 0.1. (b) Comparison for the struc- 
ture factors (lines: theory; points: simulation) for 
the same system at the thermodynamic points of (a). 



B. Systems displaying reentrant melting 

We now turn our attention to the opposite case, namely 
pair potentials belonging to the Q + -class. As an example 
within the FDM family, we have taken the model with 
parameter £ = 0.6 and performed a comparison between 
MC and MFA results. A characteristic example is shown 
in Fig. ^. As can be seen in Fig. fl](a), unlike the case 
of (5 ± -class potentials, the radial distribution function is 
completely deprived of any structure, although the ther- 
modynamic parameters are in the same regime as those 
presented in Figs. ||, || and [| In fact, in the present case, 
g(r) has a minimum at r = 0, not a maximum. This com- 
plete lack of structure is reflected in the shape of S(Q), 
shown in Fig. |^(b). 

These characteristic features for the Q + -class are not 
an artifact of the relatively high temperature chosen in 
the results of Fig. ^. They persist even at extremely 
low temperatures, provided the density is high enough. 
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FIG. 6: The quantity Airr^gir) within the diameter a of 
the PSM along the freezing line t = r). All the curves 
converge to a single one at high densities, indicating that 
the integral Airp f° r 2 g(r)dr scales linearly with density. 



This has been amply demonstrated recently for the case 
of the Gaussian core model, another member of the Q + - 
class |22| . In order to stress this point we present in Fig. 
| the g(r) and S(Q) of the GCM at t = 0.01 and r) = 6.0. 
Though g(r) displays some structure up to r w 2er, the 
structure factor S(Q) shows no signature of some kind of 
ordering [ pZi| . At any arbitrarily small but finite temper- 
ature, a high enough density can be found for which the 
MFA is valid and then the assumption that a uniform 
phase exists leads consistently to a fluid which has ideal- 
gas behavior, i.e., vanishingly small correlations. These 
liquids are different from usual ideal gases in that, e.g., 
their pressure P and isothermal compressiblity \t scale 
respectively as P ~ p 2 and \t ~ tp~ 2 . Nevertheless, 
they are thermodynamically stable. Hence, for potentials 
in the Q + -class, the equilibrium phase for sufficiently 
high densities at arbitrarily small but finite temperatures 
is the uniform fluid. 



V. GENERIC PHASE DIAGRAMS 

We now turn to the opposite limit of the low 
temperature-low density part of the phase diagram. 
There, following the original ideas of Stillinger [Q, a HS 
mapping can be performed, as follows. The Boltzmann 
factor exp[— j3v{r)\ of the potential varies monotonically 
from the value exp(— j3e) = (since j3e 3> 1 there) at 
r = to unity at r — > oo and has a close resemblance 
to that of a hard sphere system. We can thus define an 
effective hard sphere diameter ohs through the relation: 



exp[-/3v(ans)] = 1/2. 



(14) 



Writing v(r) = ecf)(r/cr) and using the fact that <f>(x) is a 
monotonic function in order to establish that the inverse 
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FIG. 7: (a) Comparison between theory (thick lines) 
and simulation (thin lines) results for g(r) of a FDM 
system with £ = 0.6. (b) Comparison for the struc- 
ture factor S(Q) (lines: theory; points: simulation) for 
the same system at the thermodynamic point of (a). 

function (f)~ 1 {x) exists, we can rewrite eq. (|l4| ) as: 

a HS =a^ 1 (t\n2). (15) 

We now use the known fact hard spheres freeze at 77ns — 
0.5 together with eq. (|l5| ) above in oder to obtain the low 
temperature-low density freezing line of the system as: 

t i (v) = ^((^r 1/3 )- (i6) 

As the limit <f>(x) — > is attained for x — ► 00 only, it fol- 
lows that the low temperature-low density freezing line 
of the systems goes to 77 = at tf = Q. Eq. ( |16| ) is 
valid for all potentials we consider here; however, for Q + - 
potentials, combining the HS-like freezing at low temper- 
atures and densities with the fact that at high densities 
the fluid has to be stable, derived in the preceding sec- 
tion, we can draw the conclusion that such systems must 
display reentrant melting and an upper freezing temper- 
ature. 



FIG. 8: Same as Fig. but for the Gaussian core model. 



We have now taken eq. ( |l6| ) for the low-i and low-77 
freezing line of the FDM and combined it with the ana- 
lytic expression at the opposite limit, eq. (0), in order 
to draw schematically the evolution of the phase diagram 
of the FDM as a function of £, for £ < £ c . The results 
are shown in Fig. ^[ With increasing £, the slopes of the 
high-i freezing lines decrease; at the limit £ — ► 0, corre- 
sponding to the PSM, the \ow-t freezing line approaches 
the horizontal axis vertically, as is dictated by the fact 
that the PSM becomes equivalent to the HS system there 
H . In the inset of Fig. ||, we show the phase diagram for 
a system with £ = 0.6 > £ c , showing reentrant melting 
behavior. The evolution of the phase diagram from a 
clustering to a reentrant melting behavior can be easily 
visualized from this picture. 

Finally, it is important to point out that Stillinger has 
proven that any system interacting by means of a po- 
tential which (i) is differentiable at least four times, (ii) 
vanishes strongly enough at infinity to be integrable and 
(iii) is +1 at the origin, will inevitably lead to a reen- 
trant melting phase diagram under the assumption that 
the competing crystal structures have single lattice site 
occupancy M. We can therefore now complete the state- 
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FIG. 9: The evolution of the phase diagram of Q ± -FDM's 
with £. To the right of the freezing lines the system is 
solid and to the left fluid. Inset: the phase diagram of 
a Q+-FDM with £ > £ c , obtained by solving the HNC 
and employing the Hansen- Verlet criterion. Below the 
bell-shaped curve the system is solid and above fluid. 

ment and say that if a potential belongs to the Q + -class, 
then it will freeze into crystals of single occupancy and 
then remelt upon increase of the density. But if it be- 
longs to the Q ± -class, then it will freeze into a clustered 
solid at any temperature. Clustering appears therefore 
to be the crucial mechanism for crystal stabilization in 
these systems. 



VI. SUMMARY AND CONCLUDING 
REMARKS 



To summarize, we have established a criterion for the 
topology of the phase diagrams resulting from repulsive, 
bounded interactions, which is very simple in its formu- 
lation and states that if the Fourier transform of the 
pair potential is positive-definite, then the system shows 
reentrant melting but if not then it freezes at all tem- 
peratures, into clustered crystals with multiply occupied 
sites. We have also established that at temperatures ex- 
ceeding the interaction strength the mean-field theory is 
reliable at all densities and its accuracy improves quickly 
with increasing temperature. We close with the remark 
that there is a certain similarity between the ideas put 
forward here and the considerations on freezing for sys- 
tems featuring of diverging interactions at infinite spatial 
dimensions J3(| [S7| |3£| [39]. However, in the latter case, 
the direct correlation function is given by the Mayer func- 
tion f(r) = exp[— (3v(r)\ — 1 of the interaction potential 
and not by —f3v{r) as in the case at hand. 
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